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Abstract. We consider a synthetic aperture radar (SAR) system that uses 
ultra-narrowband continuous waveforms (CW) as an illumination source. Such 
a system has many practical advantages, such as the use of relatively simple, 
low-cost and low-power transmitters, and in some cases, using the transmitters 
of opportunity, such as TV, radio stations. Additionally, ultra-narrowband CW 
signals arc suitable for motion estimation due to their ability to acquire high 
resolution Doppler information. 

In this paper, we present a novel synthetic aperture imaging method for 
moving targets using a bi-static SAR system transmitting ultra-narrowband 
continuous waveforms. Our method exploits the high Doppler resolution provided 
by ultra-narrowband CW signals to image both the scene reflectivity and to 
determine the velocity of multiple moving targets. Starting from the first 
principle, we develop a novel forward model based on the temporal Doppler 
induced by the movement of antennas and moving targets. We form the 
reflectivity image of the scene and estimate the motion parameters using a 
filtcrcd-backprojection technique combined with a contrast optimization method. 
Analysis of the point spread function of our image formation method shows 
that reflectivity images are focused when the motion parameters are estimated 
correctly. We present analysis of the velocity resolution and the resolution of 
reconstructed reflectivity images. We analyze the error between the correct and 
reconstructed position of targets due to errors in velocity estimation. Extensive 
numerical simulations demonstrate the performance of our method and validate 
the theoretical results. 



1. Introduction 

1.1. Motivations 

Conventional synthetic aperture radar (SAR) is designed for stationary target imaging 
[1,2]. Moving targets are typically smeared or defocused in SAR images [3]. Many 
different approaches have been suggested to address the moving target imaging 
problem for conventional SAR systems [4-22]. Both the imaging of static scenes and 
moving targets in conventional SAR rely on the high range resolution provided by 
wideband transmitted waveforms. Such waveforms are ideal in localizing the targets, 
but poor in determining their motion parameters. 

In this paper, we consider a SAR system that uses ultra-narrowband continuous 
waveforms (CW) as an illumination source. Unlike the high range resolution 
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waveforms used by conventional SAR systems, ultra-narrowband CW signals have 
high Dopplcr resolution which can be used to determine the velocity of moving 
targets with high resolution. CW radar systems also have the advantage of using 
relatively simple and low cost transmitters and receivers which can be made small 
and lightweight [23-26] . Additionally, a SAR system that uses ultra-narrowband CW 
signals may not need a dedicated transmitter. Ambient radio frequency signals, such 
as those provided by radio and television stations, etc., can be used as illumination 
sources. In [27], wc presented a synthetic aperture imaging method of stationary 
scenes using ultra-narrowband CW signals. (See also the introduction of [27] for a 
survey of stationary target imaging using Doppler only measurements.) In this paper, 
we present a new and novel method for synthetic aperture imaging of both stationary 
scenes and multiple moving targets using such waveforms. Our method exploits the 
high Doppler resolution provided by such waveforms to form high resolution images 
of both the stationary scatters and to determine the velocity of moving targets. To 
the best of our knowledge, our method is the first in the literature that addresses 
the synthetic aperture imaging of moving targets using ultra-narrowband continuous 
waveforms. 

1.2. Related Work 

Conventional wideband SAR moving target imaging techniques can be roughly 
categorized into two classes depending on the assumptions made on the motion 
parameters. 

The first class of techniques either assume a priori knowledge of target motion 
parameters or estimate this information prior to image reconstruction [4-11]. These 
motion parameters, which include relative velocity of targets with respect to antennas, 
Doppler shift, Doppler rate etc. are then used to reconstruct "focused" reflectivity 
images. However, in practice a priori knowledge of motion parameters is either 
unavailable or difficult to determine. Therefore a great deal of effort has been devoted 
to develop techniques that do not require a priori knowledge of unknown motion 
parameters for image formation [12-21]. In this class of techniques, the estimation of 
motion parameters and image formation process are performed jointly. Our approach 
falls into this class of methods where we couple the estimation of multiple target 
velocities with the reconstruction of scene reflectivity. 

Techniques for the joint estimation of motion parameters and SAR image 
formation are based on variety of approaches. These include adaptation of inverse 
synthetic aperture imaging type methods [12], [17], [28]; autofocus type methods [18], 
[16]; the keystone transform [13]; time- frequency transform based imaging methods 
[14,15]; and generalized likelihood ratio type of ideas where the reflectivity images 
are formed for a range of hypothesized motion parameters from which the unknown 
motion parameters are estimated while simultaneously forming focused reflectivity 
images [18-21]. 

1.3. Overview and Advantages of Our Work 

We note that all the work in SAR imaging of moving targets has been developed 
for the conventional wideband SAR [4-22]. Our work differs significantly from the 
existing work in SAR imaging of moving targets. Conventional SAR moving target 
imaging methods ignore the "temporal" Doppler since wideband waveforms have poor 
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Dopplcr resolution. Furthermore, they rely on start-stop approximation [1,2]. We 
instead begin with the wave equation and derive a novel forward model that includes 
temporal Dopplcr parameters induced by the movement of the antennas and moving 
targets. Next, we develop a novel filtered-backprojection type method combined with 
image contrast optimization to reconstruct the scene reflectivity and to determine the 
velocity of moving targets. 

Similar to [18-21], we adopt a generalized likelihood ratio type approach and form 
a set of reflectivity images for a range of hypothesized velocities for each scatterer. 
Our imaging method exploits high Doppler resolution of the transmitted waveforms 
in that we form reflectivity images by filtering and backprojecting the preprocessed 
received signal onto the position-space iso-Doppler contours defined in this paper. The 
scatterers that lie on the position-space iso-Doppler contours can be determined with 
high resolution due to high resolution Dopplcr measurements. We show that when the 
hypothesized velocity is equal to the correct velocity of a scatterer at a given location, 
the singularities of the scene are reconstructed at the correct location and orientation. 
We design the filter so that the the singularities of the scene are reconstructed at 
the correct strength whenever the hypothesized velocity is equal to the true velocity 
of a scatterer. This filter depends not only the antenna beam patterns, geometric 
spreading factors etc., but also the hypothesized target velocity. We next use the 
contrast of the reflectivity images to determine the velocity of moving targets. We 
present the point spread function (PSF) analysis and the resolution analysis of our 
method. The PSF analysis shows that our reflectivity image reconstruction method 
uses temporal Doppler and Doppler-rate in forming a high resolution image. We 
analyze the resolution of the reconstructed reflectivity images and the resolution 
of achievable velocity estimation. Our analysis identifies several factors related to 
the imaging geometry and the transmitted waveforms that effect the resolution of 
reflectivity images and velocity field. We analyze the error between the correct and 
reconstructed positions of the scatterers due to error in the hypothesized velocity. 
We derive an analytic formula that predicts the positioning errors/smearing caused 
by moving targets in reflectivity images reconstructed under the stationary scene 
assumption. Specifically, we show that small errors in the velocity estimation results in 
small positioning errors in the reconstructed reflectivity images. We present extensive 
numerical simulations to demonstrate the performance of our method and to validate 
the theoretical findings. 

In addition to the advantages provided by the ultra-narrowband CW signals, 
our moving target imaging method also has the following advantages as compared 
to the existing SAR moving target imaging methods: (1) Unlike [7-13,16-20], our 
method can reconstruct the images of multiple moving targets regardless of the 
target speed, the direction of target velocity and target location; and determine the 
two-dimensional velocity of ground moving targets. Furthermore, our method can 
reconstruct high-resolution images of stationary and moving targets simultaneously. 
(2) Unlike [4-11], our imaging method does not require a priori knowledge of the 
target motion parameters. Furthermore it does not require a priori knowledge of the 
number of moving targets present in the scene. (3) Our method focuses moving targets 
at the correct locations in the reconstructed reflectivity images. The localization and 
repositioning techniques of moving targets used in most conventional SAR or ground 
moving target indicator methods are not needed [16,18-20]. (4) Wc use a linear model 
for the target motion. However, our method can be easily extended to accommodate 
arbitrary target motions, such as nonlinear, accelerating targets. (5) It can be used 
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for arbitrary imaging geometries including arbitrary flight trajectories and non-flat 
topography. Furthermore, our image formation method is analytic which can be 
implemented computationally efficiently [29-31]. 

1-4- Organization of the Paper 

The remainder of the paper is organized as follows: In Section 2, we present our 
moving target, incident and scattered field models, and the received signal model 
from a moving scene. In Section 3, we develop a novel forward model that maps 
the reflectivity and velocity field of a moving scene to a correlated received signal. 
In Section 4, we develop an FBP-typc image formation method to reconstruct the 
reflectivity of the scene and a contrast- maximization based velocity estimation method. 
In Section 5 we analyze the resolution of the reconstructed reflectivity images and the 
velocity resolution. In Section 6 we present the error in the position of scatterers due 
to error in the hypothesized velocity. In Section 7, we present numerical simulations. 
Section 8 concludes our paper. 

2. Model for Moving Targets, Incident Field, Scattered Field and 
Received Signal Models 

We use the following notational conventions throughout the paper. The bold Roman, 
bold italic and Roman lower-case letters are used to denote variables in R 3 , R 2 and 
R, respectively, i.e., x = (x,x) € R 3 , with aid 2 and The calligraphic letters 

(J-,IC etc.) are used to denote operators. Table I lists the notations used throughout 
the paper. 

Let the earth's surface be denoted by x = (x,ip(x)) £ R 3 , where x € R 2 and 
tjj : R 2 —¥ R is a known function for the ground topography. Furthermore, we assume 
that the scattering takes place in a thin region near the surface. Thus, the reflectivity 
function has the form 



2.1. Model for a Moving Target 

Let z = r(x, t) denote the location of a moving target at time t, where x denotes the 
location of the target at some reference time, say t = 0. We assume that for each 
t E [0,T], the function T(-,t) : R 3 — >• R 3 is a diffeomorphism. Physically, this means 
that two distinct scatterers cannot move into the same location. Furthermore, we 
assume that for each x £ R 3 , T(x, •) : R — ► R is differentiable. 

Let the inverse r _1 (-,£), of the function T(-,t) be a(-,i), i.e., x = a.(z,t) = 
r _1 (z,i). We assume that the refractive indices of the scatterers are preserved over 
time, however, the scattercr at x moves along the trajectory z = T(x, t). Thus, V(x) 
at time t = translates as V(a(z,t)) at time t. 

Let v x (i) denote the velocity of the target at time t, located at x when t = 0, i.e., 



where T(x, t) = [Ti(x,t), ^(x, t), T3(x, t)] T and f(x, •) denotes the derivative of 
r(x, •) with respect to t. We define 



^(x) = p(x)S{x - ^i(sc)). 



(1) 



v x (t) = r(x,t) 



(2) 



[r 1 (x,t),r 2 (x,t),r 3 (x,t)] 



v x (t) = [Ti(x,*), r 2 (x,t)]. 



(3) 
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Table 1: Table of Notations 



Symbol 


Designation 


Wo(/o) 


(Angular) carrier frequency of the ultra-narrowband waveform 


x = (x,ip(x)) 


Earth's surface 


V(x) 


3D Reflectivity function 


p{x) 


Surface reflectivity 


r(x,t) 


Location of the moving target at time t located at x at t = 


V x 


Velocity of the moving target located at x at time t = 


X 


Unit vector in the direction of x 6 R 3 




Flight trajectory and velocity of the transmitter 


7h(*)»7h(*) 


Flight trajectory and velocity of the receiver 


r(t) 


Received signal along the receiver trajectory ~fR(t) due to a 




transmitter traversing the trajectory 7r(i) 


g(as,«) 


Reflectivity function of the moving target that takes into account 




the target movement 


p(*),p(*) 


Transmitted waveform and its complex amplitude 


s 


Temporal translation variable 


A 4 


Temporal scaling factor 




Temporal windowing function 


d(s,/i) 


Windowed, scaled-and-translated correlations of the received signal 




and the transmitted waveform 


J" 


Forward modeling operator 


ip{t,x,v,s,ii) 


Phase of the operator T 


A(t, x, v, s, n) 


Amplitude of the operator T 


supp(A) 


Support of A 


/ d (s,x,v) 


Bistatic Doppler frequency with respect to a moving target 


F(a,At) 


Four-dimensional bistatic iso-Doppler manifold 


F„ (s,/x) 


Two-dimensional position-space bistatic iso-Doppler contours 


F Xo (s,/x) 


Two-dimensional velocity-space bistatic iso-Doppler contours 


/dO,x, v) 


Bistatic Doppler-rate with respect to a moving target 


F(s,C) 


Four-dimensional Bistatic iso-Doppler-rate manifold 


F, (s,C) 


Two-dimensional position-space bistatic iso-Doppler-rate contours 


Fx (s,C) 


Two-dimensional velocity-space bistatic iso-Doppler-rate contours 




b utered-backprojection reflectivity imaging operator tor a hypoth- 




esized velocity v^. 


LZl(z,x) 


Point spread function of K, Vh 


Q Vh (z,t,s) 


Reconstruction filter of the reflectivity imaging operator IC Vfl 


Pv h (z) 


Reconstructed reflectivity image for a hypothesized velocity 




Data collection manifold at x = (x, ip(x)) for v; t = vo 




Length of the support of the temporal windowing function <j>(t) 


I{Vh) 


Contrast-image 


M 


Sample mean over the spatial coordinates 


<; 


Fourier vector associated with the velocity 
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For ground moving targets, since 

r 3 (x,f) = v((ri(x,t), r a (x,t))), (4) 

we write 

r 3 (x,t) = v„^(x) • tfi(x,t), r 2 (x,t)] 

= V.^ss) • w x (t) (5) 
where V x ip(x) is the gradient of ip{x) with respect to x. Thus, 

v x (t) = [«*(*), V x ?A(a;) • • (6) 

,2. 12. Model for the Incident Field 

For a transmitter with isotropic antenna located at z transmitting a waveform s(t), 
the propagation of electromagnetic waves in a medium can be described using the 
scalar wave equation [32,33], 

[V 2 -±d?]E(t,y)=5(y-z)s(t) (7) 
cr 

where c is the speed of electromagnetic waves in the medium and E(t, x) is the electric 
field. Note that this model can be extended to include realistic antenna models in a 
straightforward manner. 

The propagation medium is characterized by the Green's function, which satisfies 

tv a -ia?Mt >y ) = -*(y)*(t). (8) 



In free-space, the Green's function is given by 

/ A *(*-|y|/co) (Q , 

g(y ' t)= 47 r|y| (9) 

where Co is the speed of light in vacuum. 

Let -frit) be the trajectory of the transmitter and p(t) be the transmitted 
waveform. The incident field E m satisfies the scalar wave equation in (7) where c 
is replaced by Co and z is replaced by 7r(^) : 

[V 2 - \d?]E™(t, z) = 6(z - lT {t)) P {t). (10) 
c o 

Thus, using (9), we have 
2. 3. Models for the Scattered Field and the Received Signal 

Let E sc (y,t) denote the scattered field at y due to the transmitter located at -frit) 
transmitting waveform p(t) . Then, using (7) and under the Born approximation and 
the assumption of isotropic receiving antenna, we have 
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Let 7_r(£) denote the trajectory of the receiver and r(t) denote the received signal 
at the receiver. Then, we have 



r(t) = E sc ( lR (t),t) 

S(t -f- \ lR (t) - z|/c ) 

47r|7^(t) - z| 
S(t'-t"-\z- lT (t")\/c ) 



V( a (z,t')) 

p(t")dt"dt'dz. (13) 
4tt\z — ■jryt )\ 

Assuming that the waveform is transmitted starting at time s, for a short duration 
of t" G [0,T] t, the wave goes out at s + 1", t" G [0,T] from the transmitter, reaches 
the target at t' + s and arrives at the receiving antenna at t + s. Note that t", t', t arc 
relative time variables within the interval that starts at time s. Thus, for this short 
time interval, using (13), we have 

f S (t-e | 7 ( . + .) zi/co) M , + 

J 47T|7fi(t + s) - z| 

x .(f-f-i-y f + „ dt , dz 

47r|z-7 T (t" + s)| 
In (14), we make the following change of variables 

z -> x = a(z,t' + s) = r _1 (z,t' + s) (15) 

and obtain 

/ - < - y ' - r y + '"^ v(x)|v,r ( x, r + .)), 

J 47r|7 fl (t + s) - T(x,t' + s)\ 

47r|r(x,t' + s)-7 T (t" + s)| ^7 V ; 

where | V x r(x, t' + s)) \ is the determinant of the Jacobian that comes from the change 
of variables. 

We make the assumption that the scatterers are moving linearly and therefore 

r(x,i)wx + v x i (17) 

where the velocity v x is now time independent. Furthermore, we assume that 
|V x r(x, t)\ « 1 since radar scenes are not very compressible. Thus, (16) becomes 

r( , +s)= rKt-t'- hR (t +S )-( X+ ,(t' + sW c ) 

^ ' J 47r| 7fi (< + .s)-(x + v x (f' + S ))| 1 ' 

x ^-^-\ X + ^ + S )-7T^ + s)\/ C0 ) ^ 

47r|x + v x (t' + s) - lT (t" + s)\ Fy ' y ' 

Note that in (18) x = [x,ip(x)] and v x = [v x , V x ip(x) ■ v x ]. 
We now define 

q{x,v) =p(x)5(v-Vx) (19) 
«p(*Mt;,w x ) (20) 

J For a typical wideband chirp pulse, this time interval is in the order of 10~ 6 s, while for an 
ultranarrowband CW signal, it is in the order of 10 — 3 s or longer. 
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as the phase-space reflectivity function of the moving scene where (p is a diffcrcntiablc 
function of v that approximates S(v — t> x ) in the limit. Using (19) and (1), we rewrite 
(18) as follows: 

1 ' J 4-ir\ lR (t + s)-(x + v(t> + s ))\ qy ' 

x 5{t ' - r = |x ^ ± s) - ± : )i/c ° + ^dxd. (2D 

47r|x + v(t' + s) -7 T (*" + s)| v ; v ; 

where v = [v, \7 x tp(x) ■ v]. 

We next make some approximations to evaluate t',t" integrals in (21). First, we 
make the Taylor series expansions in 7^ and 7t around t,t" = 0, 

lR(t + s) »7sW+7# + - 1 (22) 
lT (t" + s) « 7r(s) + 7t(s)*" + • • • . (23) 
Next, under the assumptions that 

|x + v(t' + s) - 7T (s)| » |vt'|, |7r(*)t"| 

| 7i? ( s )-(x + v(i' + S ))|>K|, | 7il ( S )i|, (24) 

we approximate 

\x + v(t' + s) - ~/ T (t" + s)\ « |x + vs-7 T (s)| 

+ (x + wP 7t (s) ■ [vi' - 7T («)t"] , (25) 
|7«(* + s) - (x + v(t' + s))| « | 7R (s) - (x + vs)| 

+ 7r(s)^T+vs) • [ 7fi (s)t - vf'] . (26) 

Thus, substituting (25) and (26) into (21) and carrying out t" and t' integrations, 
we obtain 

p(at — t + s)q(x 7 v) 
(47r) 2 | 7i? (s) - (x + vs)||(x + vs) - 7 t(s)| 
where the time dilation a is given by 

a = 1 - Jr(s) - (x + vs) • 7fl(g)/co _ l+7r(g)-(x + vs)-v/c 

1 + 7t(s) - (x + vs) • 7t(s)/c 1 - 7i?( s ) - ( x + vs ) ■ v / c o 
and the time delay r is given by 

r re [|7t(s) - (x + vs)| + |7b(s) - (x + vs)|]/c 

- [(7 T (s) :r 7x + vs)+7 i? (s)'^(x + vs)) ■ vs]/co. (29) 

We see that the time dilation term a in (28) is the product of two terms. The 
first term is the Doppler scale factor due to the movement of the transmitting and 
receiving antennas. The second term is the Doppler scale factor due to the movement 
of targets. Similarly, the delay term r in (29) is composed of two terms. The first 
term represents the bistatic range for a target located at x + vs, while the second term 
describes the range variation due to the movement of targets. 

Note that conventional wideband SAR image formation methods assume that the 
radar scene is stationary. Therefore, the Doppler scale factor due to the movement 
of targets is ignored and set to 1. Furthermore, these methods rely on the "start- 
stop" approximation [1, 2] where the movement of the antennas within each pulse 
propagation is neglected. Therefore, the Doppler scale factor induced by the movement 



r(t + s)= I j—r-^L. L _7_fZ^ ^ _ _^ dx dv (27 ) 
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of antennas is also ignored and set to 1. As a result, wideband SAR imaging methods, 
including the ones developed for moving target imaging [4-22] , ignore the time dilation 
term a in (27) and set it equal to 1 since wideband signals cannot provide high 
resolution Dopplcr measurements. 



2-4- Received Signal Model 

For a narrowband waveform, we have 

p{t) = e^pit) (30) 

where loq denotes the carrier frequency and p(t) is the complex envelope of p, which 
is slow varying as a function of t as compared to e 1JJot . 
Substituting (30) into (27), we obtain 

2 [ p(at-T + s)^ at - r+ *')q{x,v) j 
r(t + s) = -0Jq / — dxdv (31) 

where a and r are as in (28) and (29) and Gtr is the product of the geometrical 
spreading factors given by 

G T r(s,x,v) = \ jr(s) - (x + vs)||x + vs - 7t(s)| • (32) 

Note that p is a slow-varying function of time. Therefore, we approximate 
p{at) m p(t) in the rest of our discussion. Furthermore, since the speed of the 
antennas and the scatterers are much less than the speed of light, we approximate 
(28) as a w 1 + (3 where 

(3 = [ 7x (s)^T+ vs) • (v - 7x (s)) + 7ii (s)^r+ vs) • (v - y R (a))]/co . (33) 

Note that /o/3, where /o = wo/27r, represents the total Doppler frequency induced by 
the relative radial motion of the antennas and the target. We refer to — foft as the 
bistatic Doppler frequency for moving targets and denote it with /^(s, x, v), i.e., 

fd(s,x,v) = ^[ 7 t(s)'^~(x + vs) • (7r(s) -v)+7h(s)^"(x + vs) • ( 7 h(s) - v)].(34) 

CO 

Note that in (32) and (34), x = [x,ip(x)} and v = [v,V x if)(x) ■ v]. 



3. Forward Model for Moving Target Imaging 

In this section, we derive a forward model by correlating the windowed and translated 
received signal with the scaled or frequency-shifted transmitted waveform, which is a 
mapping from the four-dimensional position and velocity space to the data space that 
depends on two variables, translation and scaling factor. We use the forward model 
to reconstruct the moving targets in two-dimensional position space and to estimate 
their two-dimensional velocities. 

We define the correlation of the received signal given in (31) with a scaled or 
frequency-shifted version of the transmitted signal over a finite time window as follows: 

d{s,n) = J r{t + s)p*(p£)(f>(t)dt (35) 

for some s£l and fi € R + , where 4>(t), t e [0, TA is a smooth windowing function 
with a finite support. 
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Substituting (31) into (35), we obtain 

e itJo(a — M)*gi"o(s--r) 



d(s,p) = I /A _^ 2n — t — ZTZX UJ oP(' t ~ T + s)p*(t)q(x,v)dxdvdt. (36) 



Note that since p(t) is a slow- varying function of t, we use p(fd) w p(t) in (36). 
We define the forward modeling operator, J 7 , as follows: 

d(s,n) « T[q](s,fi) 



J e -m^^)A(t,x,v,s,fj,)q(x,v)dxdvdt (37) 



where 



<p(t, x, v, s, pi) = 27r/ t[(/i - 1) + fd(s, x, v)/f ] , (38) 

, pit - t + s)p* (tW^-^ui , N 

At, x, «, a, M - ^ , 2 /f V r °- 39 

(47r)^G T fl(s,a;,v) 

and the bistatic Doppler frequency fd is as defined in (34). 
We assume that for some tha, A satisfies the inequality 

sup \d?d^d%°dl\dZdlldllA{t,x,v,s,ii)\ <C A (l+t 2 Y m --^/ 2 (40) 

(t,/i,s,£c,u)(EW 

where U is any compact subset of R x R + xlxl 2 xl 2 , and the constant Ca depends 
on U, o;t )M , /3 S , e! 2 , £i,2- This assumption is needed in order to make various stationary 
phase calculations hold. 



3.1. Leading Order Contributions of the Forward Model 

Under the assumption (40), (37) defines J 7 as a Fourier integral operator whose 
leading-order contribution comes from the intersection of the illuminated ground 
topography, the velocity field whose third component lies on the tangent plane of 
the ground topography and (x, v) £ K 3 x M. 3 that have the same bistatic Doppler 
frequency 

We denote the four-dimensional manifold formed by this intersection as 

F(s, fx) = {(x, v) : f d (s, x, v) = / (1 - /x), (x, v) £ supp(A)} (41) 

and refer to F(s, p) as the bistatic iso-Doppler manifold. 

In order to visualize the four-dimensional bistatic iso-Doppler manifold for moving 
targets, we consider the cross-sections of the bistatic iso-Doppler manifold for a 
constant velocity and a constant position. We define 

F v (s,(J.) = {x : / d (s,x,v ) = /o(l - p), (x,v ) G supp(A)} (42) 

and 

F Xo (s,fj) = {v : / d (s,x ,v) = / (1 - n), (x ,v) £ supp(A)} . 

(43) 

Fig. 1 and Fig. 2 show the position-space and velocity-space bistatic iso- 
Doppler contours for three different flight trajectories over a flat topography: (a) 
The transmitter and receiver are both traversing straight linear flight trajectories. 
7t(s) = [3.5, s, 6.5]km and 7.r(s) = [( s — 7),s,6.5]km where s = vt with speed 
v = 261 m/s. (b) The transmitter is traversing a straight linear flight trajectory, 
1r(s) = [s,0, 6.5]km and the receiver is traversing a parabolic flight trajectory, 
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Figure 1: Position-space bistatic iso-Dopplcr contours determined for a certain s 
and a fixed v = [—150, 150, 0]m/s for three different transmitter and receiver flight 
trajectories indicated by the dashed and dash-dot lines, respectively. The black and 
white triangles denote the corresponding positions of the transmitter and receiver. 
Note that each red curve corresponds to a distinct value of fi. 




(a) (b) (c) 



Figure 2: Velocity-space bistatic iso-Doppler contours determined for a certain s and 
a fixed Xo = [5, 10, 0]km for three different transmitter and receiver flight trajectories 
shown in Fig. 1. Note that each red curve corresponds to a distinct value of fi. 



7T (s) = [s, (s - ll) 2 * 22/121, 6.5]km where s = vt with speed v = 261 m/s. (c) The 
transmitter and receiver are both traversing a circular flight trajectory. 7t(s) = 7c (s) 
and 7_r(s) = Jc(s— 7r/4) where 7c(s) = [11 + 11 coss, 11 + 11 sins, 6.5]km with s = j^t 
where speed v — 261 m/s and radius R = 11km. 

4. Image Formation 

A natural choice to form phase-space reflectivity images would be to use a filtered-back 
projection (FBP) type imaging operator that filters and backprojects the data onto the 
four-dimensional bistatic iso-Doppler manifolds introduced in Section 3. Ideally, we 
wish to reconstruct a phase-space reflectivity image so that the point spread function 
of the imaging operator is an approximate Dirac-delta function in both position and 
velocity spaces. However, since the data is two-dimensional and the phase-space 
reflectivity is four-dimensional, it may not be possible to obtain such a point spread 
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function by backprojccting onto the four-dimensional bistatic iso-Dopplcr manifolds. 

Therefore, we assume that the velocity is constant, say v^, and reconstruct 
a set of two-dimensional reflectivity images in position space only for a range of 
hypothesized velocities. We refer to each image as the -reflectivity image and 
form it by an FBP-type imaging operator, where we filter and backproject the data 
onto the position-space bistatic iso-Dopplcr contours, i.e., the cross sections of the 
bistatic iso-Dopplcr mainfolds for a range of hypothesized velocities. We show that 
whenever the hypothesized velocity is equal to the correct velocity for a scatterer, 
the scatterer can be reconstructed at the correct location in the position space. We 
design the FBP filter to ensure that the reconstructed reflectivity for a scatterer has 
the correct strength whenever the hypothesized velocity is equal to the true velocity 
of the scatterer. From this set of images, we estimate the velocity of the scattercrs 
using a figure of merit that measures the degree to which the images are focused. The 
reflectivity images corresponding to the estimated velocities provide focused images 
of the moving scatterers present in the scene. 

Below we introduce the FBP operator in forming the u^-reflectivity images, 
analyze its point spread function, and next present the design of the FBP filter. 
Finally, we describe how to determine the velocity of moving targets. 

4-1- Vh~ Reflectivity Image Formation 

We form the i^-rcflcctivity image q Vh {z) for a fixed hypothesized velocity = 
[vh, ^ z4>{z.) ■ Vh] by filtering and backprojecting the data onto the position-space 
iso-Doppler contour F Vh (s,fi): 

q Vh (z) := K, Vh [d\(z) 

= J c i ^^ z ' s ^ ) Q Vh {z,t 7 s)d(s,fi)dtdsdfi, (44) 

where K. Vh is the filtered-backprojection operator for the fixed velocity v^, 

<f> Vh (t,z,s,n)=<t>(t,z,v h ,s,n) (45) 

and Q Vh is the filter to be determined below. Note that Vh is a fixed parameter for 
<p Vh and Q Vh . 

We assume that for some fn,Q^, h , Qv h satisfies the inequality 

sup \d?dP>dl\dllQ Vh {z,t, S )\ <C Q ^(l+t 2 ) (m ^- M)/2 (46) 
(t,s,z)eu 

where U is any compact subset of R x K + x M 2 , and the constant Cq depends 
on IA, at,/3 s , 61,2. Under the assumption (46), (44) defines K. Vh as a Fourier integral 
operator. 

4.2. PSF Analysis 

Substituting (37) into (44), we rewrite (44) as 
q Vh (z) := lC Vh F[q](z,v h ) 

= J LZl(z,x)q(x,v*)dx (47) 

where L%?(z, x) is the Point Spread Function (PSF) of the two-dimensional reflectivity 
imaging operator for the hypothesized velocity Vh with respect to the true velocity v x 
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given by 

L Vx (z x) = J e'^'i ( t < z < s <f l )-<t>( t < 3! < v x< s <i 1 )] 

x Q Vh (z,t,s)A(t',x,v x ,s,fi)dtdsdfidt'. (48) 

We define 

*fc = <t> Vh (t, z, s, n) - <f>(t, x, v x , s, fl) 

= 27rt[(M - l)/o + fd(s, z, v h )} - 2-Kt'[{n - l)/ + f d (s, x, v x )} . (49) 

Applying the stationary phase theorem to approximate the t' and \i integrations 
in (48) f, we obtain 

d t ^ k = -2n[(p - l)/ + f d (s,x,v x )} = 
, fd(s,x,v x ) 

=^> /i = 1 , (50) 

Jo 

d^ k = 2ir(t - t')f = 
=>t = t' . (51) 
Substituting the results back into (48), we get the kernel of the image fidelity operator 



L%(z,x)k J e 



i27rt[f d (s,z,v h )-f d (s,x,v x )] 



x Q Vh (z,t, s)A(t,x,v x ,s,l - f d (s,x,v x )/ f )dtds . (52) 

To simplify our notation, we let 

A(t, x, v x , s) = A(t, x, v x , s, 1 - f d (s, x, «x)//o) • (53) 

The main contribution to L"* comes from the critical points of the phase of K, Vh J- 
that satisfy the conditions [34]: 

d t (2irt[f d (s, z, v h ) - fd(s, x, v x )]) = 

=> fd(s,z,v h ) = fd(s,x,v x ) , (54) 
d s {2nt[f d (s,x,v x ) - fd(s,z,v h )}) = 

f d (s,z,v h ) = f d (s,x,v x ) (55) 

where f d (s, x, v x ) denotes the first-order derivative of f d (s, x, v x ) with respect to time 
s, i.e., f d (s, x, v x ) = df d (s, x, v x )/ds. We refer to f d (s, x, v x ) as the bistatic Doppler- 
rate. 

Using (34), we obtain 

1 



f d {s,x,v x ) = — 

CO 



\jt(s) - (x + v x .s) 
1 



|(7t(s) -v x Ur + (7t(s) - (x + v x s)) -7 T (s) 



1(7*00 " VxUr + (7b(«) - ( x + vx*)) • Jr{s) 



(56) 



|7«(s) - (x + v x .s 
where 

Ht,r{ s ) ~ v x)± = (jTJi(s) - V x ) - 

(1t,r( s ) - ( x + v xs))[(7t,a(s) - (x + v x s)) • (7t,r(s) - v x )] (57) 

f The determinant of the Hessian of <t>t, is (2n) 2 f$. Thus, the stationary points are non-degenerate. 
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X (km) X (km) X (km) 



(a) (b) (c) 

Figure 3: Position-space bistatic iso-Doppler-rate contours determined for a certain 
s and a fixed v = [—150, 150, 0]m/s for three different transmitter and receiver 
flight trajectories as described in 3.1. The black and white triangles denote the 
corresponding positions of the transmitter and receiver. Note that each blue curve 
corresponds to a distinct value of C . 



denotes the projection of the relative velocity 7T,i?,( s ) — v x onto the plane whose 
normal vector is along Jt.r( s ) — ( x + v x s). Note that in (56) x = [x,tp(x)] and 
v x = [v x , V x ip(x) ■ v x ]. 

In (56), the summation of the first two terms in the square bracket corresponds 
to the relative radial acceleration between the transmitter and the target located at 
x + v x s at time s, while the summation of the last two terms in the square bracket 
corresponds to the relative radial acceleration between the receiver and the target 
located at x + v x s at time s. For the derivation of (56), see Appendix A. 

We refer to the locus of the points formed by the intersection of the illuminated 
surface, [x,ip{ x )], the velocity field, [v, V x ip(x) ■ v], and the set {(x, v) £ I 3 x I 3 : 
fd(s,z,v) = C}, for some constant C, as the bistatic iso-Doppler-rate manifold and 
denote it by 

F(s, C) = {{x, v) : f d (s, x, v) = C, (x, v) e supp(A)} . (58) 

We consider the cross-sections of the bistatic iso-Doppler-rate manifold for a 
constant velocity and a constant position and define 

F Vo (s, C) = {x: f d (s, x, v Q ) = C, (x, v ) & supp(A)} (59) 

and 

F Xo (s,C) = {v : / d (s,xo,«) = C, (x ,v) G supp(A)} . (60) 

(59) specifies an iso-Doppler-rate contour in the two-dimensional position space. We 
refer to this contour as the position-space bistatic iso-Doppler-rate contour for moving 
targets. Similarly, (60) specifies an iso-Doppler-rate contour in the two-dimensional 
velocity space. We refer to this contour as the velocity-space bistatic iso-Doppler-rate 
contour for moving targets. 

Fig. 3 and Fig. 4 show the position-space bistatic iso-Doppler-rate contours and 
velocity-space bistatic iso-Doppler-rate contours for three different flight trajectories 
over a fiat topography that are described in Section 3.1. 

The critical points of the phase of K, Vh T that contribute to the reflectivity image 
formation are those points that lie at the intersection of the position-space bistatic 




Figure 4: Velocity-space bistatic iso-Doppler-rate contours determined for a certain 
s and a fixed x = [5, 10, 0]km for three different transmitter and receiver flight 
trajectories as described in 3.1. Note that each blue curve corresponds to a distinct 
value of C. 



iso-Doppler contours, F Vh (s,^i) and position-space bistatic iso-Doppler-rate contours, 
Fv h {s,C). For the correct velocity, i.e., Vh = f x , this intersection contributes to the 
reconstruction of the true target t,. Note that when Vh ^ v x , the points lying at the 
aforementioned intersection may lead to the artifacts in the reconstructed reflectivity 
image. 



4-3. Determination of the FBP Filter 

We determine Q Vh so that the PSF of the two-dimensional reflectivity imaging 
operator, L^(z,x) is as close as possible to the Dirac-delta function, S(z — x) for 
Vh = v x , i.e., whenever the reflectivity at z is reconstructed for the correct target 
velocity v x . We assume that at the correct target velocity, the flight trajectory and 
the illumination pattern are chosen such that the only contribution to L£*(z, x) comes 
from those points z = x. 

Thus, we linearize /d(s, z, Vh) around z = x for Vh = « x and approximate 

f d (s,z,v h ) - f d (s,x,v h ) x>(z-x)- V z fd{s,z,v h ) . (61) 

We write 

A(t,z,v h ,s) w A(t,x,v x ,s) . (62) 
Thus, (52) becomes 

LZ h h (z,x) = J ^"- !B > B ^ s ^Q Vh (z,t,s)A{t,z,v h ,s)dtds (63) 

where 

S Vh (s,z)=2wV z f d (s,z,v h ). (64) 
For each z, we make the following change of variables: 

(t,s)^Z = t3 Vh (s,z) (65) 



| We assume that the flight trajectory and the illumination patterns are chosen such that the 
intersection has a single element avoiding any right-left type of ambiguities. 
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and write (63) as follows: 

L%(z,x)= [ j(*- mH Q Vh {z,t)A(z,v h ,t)ri(z,v h ,t)dt (66) 



where 



and 



A(*,t; h ,€) = A(t(€),*,« fc ,«(0) 



(67) 



det 



9 s S Vh (s,z) 



(69) 



is the determinant of the Jacobian that comes from the change of variables given in 
(65). 

The domain of integration in (66) is given by 



0, 



{£ = ta Vh (s, z) | A(t, z, v h ,s) ^ 0, t,se 



(70) 



We refer to Q Vh ,z as the daia collection manifold at z for Vh — v x . This set determines 
many of the properties of the reconstructed reflectivity image when Vh = v x . 
Using (64) and (34), we obtain 



where 



27t/q 

CO 



D 



D 2 = 



[D + D 2 s] 



(7t(s) - v h )± 



Hr{s) - v h )_ 



7t(s) - (z + v h s)\ hn(s) - (z + v h s) 
D 2 ■ [( 7T (s) + v h s)) + (~f R (s) v h «))]} (71) 



1 d^(z)/dZ! 

1 dip(z)/dz 2 







S^«« 



Vh2 



(72) 



(73) 



and (7T,fl — v^)j_ is the projection of 7t,_r — v), onto the plane whose normal is 
Jt,r(s) — (z + Vfts) as defined by (57). Note that v/, = [w^, V^VC 21 ) ' v h], v h = 
[vhi,Vh2\- For the derivation of (71), see Appendix B. 

To approximate the point spread function L%*(z,x) in (66) with the Dirac-delta 
function, we choose the filter as follows: 

Xa„ h ,* A*(z,v h ,£) 



Qv h (z,£) 



(74) 



where xsi„ h is a smooth cut-off function that prevents division by zero in (74). 

With this choice of filter, the resulting FBP operator can recover not only the 
correct position and orientation of a scatterer, but also the correct reflectivity at x 
whenever Vh = v x in the v^-reflectivity image. 
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4-4- Determination of the Velocity Field 

The filtered-backprojection of data results in a set of reflectivity images q Vh in the 
two-dimensional position space corresponding to a range of velocity values that is 
suitably chosen for ground moving targets. When the hypothesized velocity Vh is 
equal to the correct velocity d x , the corresponding ^^-reflectivity image is focused at 
x. We measure the degree to which the reflectivity images are focused with the image 
contrast measure [18,28] and generate a contrast-image as follows: 

,(„„) . Mv^-M^m 

\M[q Vh )Y 

where Vh = [vhi , Vh.2] is the index of the contrast-image and A4[-] denotes the sample 
mean over the spatial coordinates. Note that the image contrast can be viewed as 
the ratio of the standard deviation to the mean of the i^-rcflcctivity image. This 
figurc-of-merit was previously used in [18, 28] to determine target velocities from a 
stack of images for the conventional SAR moving target imaging. 

If there are multiple moving targets with different velocities in the scene, the 
contrast-image could have several peaks each one corresponding to the velocity of a 
different moving target. We accordingly detect the moving targets and determine their 
velocities by detecting the local maxima in the contrast-image I(vh). A threshold can 
be used in the detection, which may be determined using the Constant False Alarm 
Rate (CFAR) criterion [35]. 

In practice, the discretized and estimated velocity may deviate from the true 
velocity. In the following two sections, we analyze the velocity resolution and the 
error in the reflectivity image reconstruction due to error in the estimated velocity. 



5. Resolution Analysis 

In this section, we analyze the resolution of reconstructed reflectivity images and the 
velocity resolution available in the collected data. Our resolution analysis results are 
consistent with the Doppler ambiguity theory of ultra- narrowband CW signals [36] . 



5.1. Resolution of Reflectivity Images at the Correct Target Velocity 

To determine the resolution of the reconstructed reflectivity images, we analyze the 
bandwidth of the PSF associated with the image fidelity operator JC Vh T at the correct 
target velocity. 

Substituting (74) into (66) and the result back into (47), we obtain 
q Vh (z) := K Vh F[q](z) 

= [ S z - X] tq(x,v x )dxd£. (76) 

(76) shows that the image q Vh (z) is a band-limited version of q whose bandwidth is 
determined by the data collection manifold &v h ,z whenever the hypothesized velocity 
is equal to the true velocity. The larger the data collection manifold, the better the 
resolution of the reconstructed reflectivity image becomes. Furthermore, as indicated 
by (70) and (71), the band- width contribution of £ to the reflectivity image at z is 
given by 

(■Jt(s) - VfoU (inXs) - v/Q-L 

\jt(s) - (z + v h s)\ \jr(s) - (z + v /l s)|_ 



2tt/o 

Co 



[D 



D 2 s] 
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Table 2: Parameters that affect the resolution of the u^-reflectivity image 



Parameter 



Increase (t) 



Resolution 



Carrier frequency: fo f 

Length of the windows f 
Distance \jr(s) - (z + v h s)\, \~fr(s) - (z + v h s)| f 
Relative velocity -yr — vj, or 7_r — vj, 

Bistatic angle #t.r t 

Ground topography variations f 

Number of s samples f 



t 
t 



t 
t 

t 

t 
t 



Higher (|) or Lower (4.) 



+2 cos 



6 T R(z,v h ,s) 



D 2 ■ b TR {z,\- h , s) 



(77) 



2 



where denotes the length of the support of <p{t), 6r^(z,v^,s) denotes the unit 

vector in the direction of [(7t(s) — (z + V/js)) + (7h(s) — (z + v^s))] and 0tr(z, v/j, s) 
denotes the bistatic angle formed by the transmitter and receiver with respect to the 
target located at (z + v^s) at time s. D and D 2 are as described in (72) and (73). 

(77) shows that as the carrier frequency of the transmitted signal fo becomes 
higher, the magnitude of £ gets larger, which results in higher resolution reflectivity 
image of the moving target. Furthermore, (77) shows that the resolution depends on 
the range of the antenna to the moving target via the terms 7t(s) — (z + v/js) and 
\1r(s) — (z + v/is)|; and the relative speed between the transmitter (receiver) and 
the moving target via the terms ('Jt(s) — v^)j_ and ('Jr(s) — Vh)±- As the antennas 
move away from the target, or the relative speed decreases in certain directions, the 
magnitude of £ decreases, which results in reduced resolution. Additionally, larger 
number of processing windows, i.e., s samples, used for imaging leads to a larger 
data collection manifold, and hence better resolution. As indicated by the second 
line of (77), the resolution of the reflectivity image also depends on the bistatic angle 
Otr(%, Vft,, s). Larger the 9tr(z, "Vhi s), lower the resolution becomes. 

We emphasize again that this analysis holds only for those reconstructed 
scatterers at x whose velocity v x is equal to the hypothesized velocity Vh- 

We summarize the parameters that affect the resolution of the reconstructed 
moving target image in Table II. 

5.2. Velocity Resolution 

Our imaging method discretizes the range of the velocity and forms a reflectivity image 
corresponding to each velocity sample. Therefore, the velocity resolution depends on 
how finely the range of the velocity can be sampled, which, in turn, depends on the 
"velocity bandwidth" available in the data. We show that, for a point target located 
at a fixed position, the data can be interpreted as the bandlimited Fourier transform of 
the phase-space reflectivity function with respect to the velocity variable and analyze 
the bandwidth of the data in terms of the imaging geometry, parameters of the ultra- 
narrowband CW signals, and other data collection parameters. 

We assume that the scene consists of a moving point target located at Xq at t = 
moving with velocity v Xo , i.e., 



q(x, v)=8(x- x )q(x , v). 



(78) 
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Without loss of generality, we assume that v xa = 0. Performing Taylor series 
expansion in the phase of the forward model given by (38) around v = 0, we get 

4>(t, cc, v, s,fi) w (p(t, x, 0, s, (j) + V„0(i, x, v, s,fi)\ v=0 ■ v . (79) 
Substituting (78) and (79) into (37), we obtain 

F[q](s, fi) w /' e - i #( t > a oA*>*0 e - iV «^(*> SE o>«>aiM)l«=o 

x q(xo,v)A(t 7 x 0} v,s, fi)dvdt (80) 



where 
and 



JB , o, a, M) = M(M " l)/o + fd(s, x , 0)] (81) 



V v <l>(t,x Q , v,s,/j,)\ v=0 = 2TTtV v f d (s, x ,v)\ v=0 . (82) 

Note that fd(s, xq,0) in (81) represents the Doppler frequency induced by the 
movement of the transmitter and receiver. It does not depend on the target velocity. 
Substituting (81) and (82) into (80), we obtain 

.F[g](s,//)« J (e- i2vtv -M 8 ' Xo ' v ^°- v q{xo,v)A(t,Xo,v,s, f i)dv^ 

x e -i2irt[(/i-l)/o+/ d (s,xo,0)]^ i 

Let 

tVvfd(s,x ,v)\ v=0 = <;. (84) 

We see that ^ is the Fourier vector associated with v. Therefore, the length and 
direction of <; determine the velocity resolution available in the data, d(s,/z). 
The bandwidth contribution of <j is given as follows: 

kl = \tVvfd(s,x ,v)\ v =o\ 



D ■ 



6> Tfl (x ,v, s)~ 
2 cos o T r{x , v, s) 

(jr(s) - v)± s (ir(s) - v)j_ s 



|7t(s) - (x + vs)| |7 fl (s) - (x + vs)| 



(85) 



where L^, 9tr(x-o, v, s), brFi(xo, v, s) are as defined in (77). Note that xo = 
[ceo, V'(iKo)], v = V a ; V'(3 ; o) • u] and D is given by (72) with z replaced with xq. 

Comparing (85) with (77), we see that similar to the reflectivity image formation, 
the larger the carrier frequency /o and the support of <fi(t), the higher the velocity 
resolution is. Furthermore, the velocity resolution also depends on the range of the 
antennas to the moving target via the terms I'Jt(s) — (xo+vs) and 7_r(s) — (xo+vs)|; 
and the relative speed between the transmitter (receiver) and the target via the terms 
('Vt(s) — v)j_ and ('Jr(s) — v)j_. The increase in the number of s samples used for 
imaging also results in a larger data collection manifold and hence better resolution. 
Additionally, the larger the bistatic angle 9tr(x-o,v, s) is, the lower the velocity 
resolution becomes. Note that the bistatic angle 0tr(xq, v, s) has a larger impact 
on the velocity resolution than on the position resolution due to the dependence on 
D instead of D 2 . 

Note that the parameters that affect the resolution of reflectivity images and 
velocity resolution identified in our analysis are consistent with the Doppler ambiguity 
theory of ultra-narrowband CW signals [36] . 
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6. Analysis of Position Error in Reflectivity Images due to Incorrect 
Velocity Field 

In the previous section, we show that the image fidelity operator K,x, h T reconstructs 
the singularities at the intersection of the bi-static iso-Doppler and iso-Doppler-rate 
manifolds defined by the following two equations: 

fd(s,z,v h )-fd(s,x,v x ) = (86) 
f d (s, z, v h ) - f d (s, x, v x ) = (87) 

where fd = d s fd- When Vh = v x , one of the solutions of (87) is z = x which 
shows that K. Vh J- reconstructs a singularity that coincides with the visible singularity 
of the scene, q(x,v x ). We shall refer to such singularities of K. Vh J-[q] as the useful 
singularities jj. If, on the other hand, Vh ^ v x , the useful singularities of fC Vh J-[q] 
no longer coincide with the visible singularities of the scene reflectivity. In this 
section, we analyze the shift in the location of the useful singularities Vh -reflectivity 
image due to errors in the hypothesized velocity field Vh- The analysis provides the 
positioning error between the correct and reconstructed targets due to error in their 
hypothesized velocities. Additionally, it shows the geometry and degree of smearing 
in the reconstructed reflectivity images due to incorrect velocity information given the 
imaging geometry. For simplicity, we assume that the ground topography is flat for 
the rest of our analysis. 

Suppose for the target located at x at t — moving with velocity v x , we use an 
erroneous hypothesized velocity 

Vh =v z + eAv z 



in the backprojection, where v z = v x and eAv z , e G R, is the error in the velocity Vh- 
Then, the target at position x = z is reconstructed at z e = z + Az and we have 

fd{s, z + Az, v z + eAv z ) - f d (s, x, v x ) = (89) 
f d (s, z + Az, v z + eAv z ) - f d (s, x, v x ) = 0. (90) 

(89) and (90) show that the visible singularity at x in the scene is mapped to a 
singularity at z e = z + Az in the reconstructed image. 

We want to determine the first order approximation to the shift Az due to the 
velocity error eAv z . In order to determine Az, we assume that e — s> is small and 
expand (89) and (90) in Taylor series around e = and keep the first-order terms in 
e. Then, using (86)- (87) and (89)- (90) in the Taylor series expansion, we obtain 

edj d (s, z, v z + eA-y z )| e=0 + V*/ d (s, z, v z ) ■ Az = (91) 
edj d {s, z, v z + eA-i; z )| e= o + V z / d (s, z, v z )-Az = 0. (92) 

Evaluating (91) and (92) for the bi-static Doppler frequency of moving targets, 
(91) simplifies to 

(7t(s) - v z ) 



-esA^ T 
-esAv^ R 



\jt(s) - (z + v z s)| 
\jr(s) - (z + v z s) 



J Note that in addition to useful singularities, K Vh T may reconstruct additional artifact singularities 
that are of the same strength as the useful singularities. The location of these singularities are given 
by the solution of (87) when = v x . 
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-eAv 
= Az ■ S Vm (s, z 
where for flat topography, 



(7t(s) - (z + v z s) + (■Jr(s) - (z + v z s) 
Co 



2tt/ 



(7t(s) 



(7ii(s) - v z )_l 



(93) 



(94) 



A- 



r -L,T,R 



\1t(s) - (z + v z s)| |7 fl (s) - (z + v z s) 
and (7t,.r(s) — v z )^ are the projections of Av z and 7t,_r(s) — v z onto the 



plane whose normal direction is (■Jt.r( s ) ~ (z + v z s). Note that in (93) and (94), 
v z = [v z , 0], z = [z,0], and v z = w x and z = x. In other words, z and v z are the 
correct position and velocity of the target in the image domain that is located at x 
moving with velocity v x in the scene. 
Similarly, (92) simplifies to 

sj T (s) 2(7 T (s)-v z ; 



-eAv, 



+esA 



-eAv, 



+esAv 



= -Az ■ 




I |7r(s) - (z + v z s)| 

|(7t(s)-v z )J 2 
7t(s) - (z + v z s)| 2 



»)l 17*00 
|(7ii(s) - v z 



|2 



v z s) 



C T (z, v z , s) 



■C_r(z,v z ,s) 



|7h(s) - (z + v z s) 



(95) 



where H^ z (s, z) = 9 s S Vz (s, 2;) is the derivative of S Vx (s,z) with respect to s for flat 
topography. For the explicit form of Ct,r{z>,v x ,s) and E„ z (s,z), and the derivation 
of (93) and (95), see Appendix C and Appendix D. Note that in (95), 7^ R is the 
projection of the acceleration, 7t,jj, of the transmitting/receiving antenna onto the 
plane whose normal direction is (7T r(s) — (z + v z s)). 

The shift Az in the useful singularity lies at the intersection of the solution of 
(93) and (95). (93) and (95) show that when the error in the velocity v z is in the order 
of e, the shift in the reconstructed useful singularities is also in the order of e, which 
means that the reconstructed reflectivity images would vary smoothly with respect to 
the change in the velocity around the correct value. 

(93) and (95) show that for a given aperture location V, the shift in position, 
Az, depends on the components of the velocity error Av z in the look directions of 
the transmitting/receiving antennas, (~/t,r(s) — (z + v z s)), and its projections onto 
the planes perpendicular to the antenna look directions. Clearly, the shift in position 
depends on the antenna flight trajectories. The V dependency of the shift explains 
the smearing observed in the final backprojected data. 

Clearly, (93) and (95) can be used to predict the positioning errors caused 
by moving targets in reflectivity images reconstructed under the stationary scene 
assumption, in which case Av z = — v x . 

Example - Monostotic SAR traversing a linear flight trajectory 

To understand the implications of (93) and (95) and to illustrate the shift 
in position, we consider a relatively simple scenario where the transmitting and 
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receiving antennas are colocated, traversing a linear trajectory forming a relatively 
short synthetic aperture. 

Let 7 = 7t = Jr denote the flight trajectory. Then (93) and (95) become 

A ± (7(g) ~ Vz) A f f S "T^ \ \ A frOO - V z)-L 

-esAv z ■ — — - eAv z ■ {■y(s) — (z + v z s)) = -Az ■ 



| 7 ( s )-( z + v zS )| " ^ ' ^ * " | 7 (s)-(z + v zS )| 

(96) 



and 



-esAv^ 
+2eAvf 



| 7 (s)-(z + v z s)| 
{j(s) - v z ) ± 



| 7 (s)-(z + v z s)| 

(7(s) - v z )_l[(7(s) - v z ) ■ 7(s) - (z + v z s)]s 



-esAv z • (7(3) - (z + v z s)) 
-Az • 7(s) - (z + v z s 

-Az • (7(5) - v. 



7(s) - (z + v z s)| 2 

l(7( S )-v z U| 2 



| 7 (s) - (z + v z s)| 2 

l(7(^-v z )±| 2 
7(3) - (z + v z s)| 2 

2(7(3) - v z ) • 7(s) - (z + v z s) 



|7(s) - (z + v z s)| 2 

+A - | 7W -^tv, 8 )| - < 97 » 

We assume that the radar flics along a linear straight trajectory with a constant 
velocity and observes a region of interest in the far-field and in the boresight direction 
of the antenna. Since the speed of the target is usually much smaller than the speed 
of the antenna, we assume that the relative velocity vector, 7 — v z , is perpendicular 
to the radar line of sight (RLOS), 7(s) — (z + v z s) throughout a short synthetic 

aperture. Under these assumptions, 7 = 0, (7(s) — v z ) • 7(s) — (z + v z s) = and 
(7 — v z )_l ~ 7 — v z . Thus, (96) and (97) reduce to 

|Az x l = e3]Av z L ]co S g + e|Av z ] |7(s) ~^ Z + VzS)l (98) 

and 

|Az r | = -es|Av^| +2e|Av z L |cos6» ,7(s) ~ (z + , VzS)l (99) 

l7-v z | 

where 9 denote the angle between Av^ and 7 — v z on the plane normal to the RLOS, 
I Az -1 - 1 = Az • 7 — v z , denotes the position shift along the direction of the vector 

7 — v z , i.c, perpendicular to the RLOS, and |Az r | = Az • (7(s) — (z + v z s)), denotes 
the position shift along the RLOS. We refer to | Az -1 - 1 and |Az r | as the tangential 
position error and the radial position error, respectively. Similarly, we refer to | Av^ | 
and |Av z | as the tangential velocity error and the radial velocity error, respectively. 

From (98) and (99), we see that for a fixed time (aperture point) s, the tangential 
position error, |Az |, mainly depends on the radial component of the velocity error, 
|Av z |, due to the range term, (7(5) — (z + v z s)|. Similarly, the radial position 
error, |Az r |, mainly depends on the tangential component of the velocity error, 
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|Av^|. Under the far- field assumption and for a short synthetic aperture, we note 

that the RLOS vector, (j(s) — (z + v z s), and the second terms in (98) and (99) are 
approximately s independent. In the following section, we elaborate on this example 
and show the shift and smearing for a point target moving perpendicular to the antenna 
flight trajectory. 

7. Numerical Simulations 

Wc performed two sets of numerical simulations to demonstrate the performance of our 
imaging method and to validate the theoretical results. In the first set of simulations, 
we numerically studied the reflectivity (or position) reconstruction performance and 
the velocity estimation performance of our method for a single point moving target. 
We also demonstrated the theoretical velocity error analysis described in Section 6 
using the experimental results of the first set of simulations. In the second set of 
simulations, we demonstrated the performance of our imaging method for multiple 
moving targets. Different transmitter and receiver trajectories were used in the two 
sets of simulations. In the first set of simulations, we considered a monostatic antenna 
traversing a straight linear trajectory. In the second set of simulations, we considered 
a bistatic setup where both the transmitter and receiver are traversing a circular 
trajectory. 

For all the numerical experiments, we assumed that a single-frequency continuous 
waveform operating at /o = ujq/2it = 800MHz is being transmitted. We used (27) and 
(35) to generate the data. We used (31) to generate the received signal and (35) to 
generate the data used for imaging and chose the windowing function cf> in (35) to be 
a Hanning function. 

7.1. Simulations for a Point Moving Target 

We considered a scene of size 256 x 256 m 2 with flat topography centered at 
[11, 11, 0]km. The scene was discretized into 128 x 128 pixels, where [0,0, 0]m and 
[256, 256, 0]m correspond to the pixels (1, 1) and (128, 128), respectively. We assumed 
that a point moving target with unit reflectivity was located at the center of the scene 
at time t = moving with velocity [0, 6.2, 0]m/s. Note that this position corresponds 
to the (65.65)th pixel in the reconstructed scene. 

We considered a monostatic antenna traversing a straight linear trajectory, 
7l(s) = (s, 0,6.5) km, at a constant speed. Hence, s = vt where v = 261m/s is 
the radar velocity. Fig. 5 shows the 2D view of the scene with the target and antenna 
trajectories. The aperture length used for the image was 5.5e3m, as indicated by the 
red line. 

We assumed that the velocity of the target is in the range of [—10,10] x 
[—10, 10]m/s and implemented the velocity estimation in two stages, each one using a 
different discretized step: We first discretized the entire velocity space into a 21 x 21 
grid with a step size of lm/s, from which we obtain an initial estimate of the target 
velocity, « Xj o. Then, we discretized a small region of size [—1, 1] x [—1, l]m/s around 
the initial velocity estimate into a 21 x 21 grid with a step size of O.lm/s to refine our 
velocity estimate obtained in the first stage. 

We reconstructed p Vh (x) images via the FBP method as described in Section 
4.1 with fo = 0.8c9GHz, = 42.67ms and the aperture sampling frequency, 
f s = 97.1869Hz. 
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Figure 5: 2D illustration of the simulation setup for a single moving target using a 
monostatic antenna. The dark region shows the scene considered. The red dot shows 
the position of the point target with the arrow indicating the direction of the target 
velocity. The radar platform traverses a straight linear trajectory, as shown by the 
black line. The aperture used for the image is shown by the red line. 



We formed the contrast images for each velocity estimation stage as described in 
Section 4.4. The results are shown in Fig. 6(a) and Fig. 6(b). The red circle shows 
the velocity estimation. The initial velocity estimate of t) x ,o = [0,7]m/s, is shown in 
Fig. 6(a), which is close to the true target velocity, v x = [0, 6.2]m/s. The bright region 
around the peak in the contrast image shown in Fig. 6(a) indicates that the image 
contrast varies smoothly with the hypothesized velocity. 

The contrast image obtained using a finer dscretization step is shown in Fig. 6(b). 
This contrast image results in a velocity estimate of [— 0.4, 6]m/s which is shown by 
the red circle. The estimate deviates slightly from the true value shown by a black 
circle. Looking at Fig. 6(b), we see that the refinement of velocity estimation is not 
as good as expected. This may be explained by the velocity resolution provided by 
the linear flight trajectory and the short aperture as well as waveform parameters. 

Fig. 7(a) shows the reconstructed reflectivity image of the moving target when 
Vh = v x = [— 0.4, 6]m/s. Note that the black circle shows the true target 
location. Fig. 7(b) shows the reconstructed reflectivity image of the target when the 
hypothesized velocity is equal to the true target velocity, i.e., Vh = v x = [0, 6.2]m/s. 
We see that the moving target in Fig. 7(a) is reconstructed almost as good as the one 
in Fig. 7(b) with the exception of slight energy spread and a position error due to 
error in the estimated velocity. In the following subsection, we present a quantitative 
numerical study to demonstrate the results of the analysis described in Section 6. 

7.2. Numerical Analysis of the Position Error due to Velocity Error 

We use the simulation results obtained in the previous subsection to demonstrate the 
theoretical analysis presented in Section 6. Note that the geometry considered in the 
simulation is consistent with the example given in subsection 6. 



Synthetic Aperture Imaging of Moving Targets using Ultra- Narrowband CW 25 




-10 -5 5 10 6 6.2 6.4 6.6 6.8 7 7.2 7.4 7.6 7.8 8 

v h2 (mis) V h2 (m/s) 



(a) (b) 

Figure 6: The contrast images formed in the two-stage velocity estimation, (a) Contrast 
image formed for the entire velocity space discretized using a step size of lm/s. The 
estimated velocity is 6 x ,o = [0, 7] m/s as indicated by the red circle, (b) Contrast image 
formed for a small region of size [—1, 1] x [—1, l]m/s around t) x , using a step size of 
O.lm/s. This corresponds to the region [—1, 1] x [6, 8]m/s. The estimated velocity at 
this stage is v x = [—0.4, 6]m/s, as shown by the red circle. The black circle shows the 
true target velocity. 

Fig. 8(a) and Fig. 9(a) show the reflectivity images reconstructed using Vh = 
[0, 6.7]m/s and Vh = [0.5, 6.2]m/s, respectively. Note that the former has a radial 
velocity error of |Av£| = 0.5m/s and the latter has a tangential velocity error of 
|Av^| = 0.5m/s. The true position of the moving target is shown by a red circle in 
Fig. 8(a) and Fig. 9(a). 

As compared to the image reconstructed using the correct velocity shown in 
Fig. 7(b), we see that in Fig. 8(a), there is an obvious horizontal (or tangential) 
position shift, while in Fig. 9(a), there is an even larger vertical (or radial) position 
shift. This is predicted by (98) and (99) in Section 6, which state that the velocity error 
in the tangential direction would lead to roughly twice the radial position error that 
would result from the same magnitude of velocity error in the radial direction. Table 
III compares the position shift errors that are measured from the reconstructed images 
and the ones predicted by (98) and (99), as well as the estimated target positions (in 
pixel indices) and the corresponding reflectivity values. We also note the smearing in 
Fig. 8(a) and Fig. 9(a) due to the velocity error. This can be seen more clearly by 
comparing the X and Y profiles of the reconstructed images, as shown in Fig. 8(b), 
Fig. 8(c) and Fig. 9(b) and Fig. 9(c). Note that the X and Y profiles were shifted to 
the center for ease of comparison with the results obtained using the correct velocity. 



7. 3. Simulations for Multiple Moving Targets 

In this subsection, we perform simulations for a scene containing multiple moving 
targets to demonstrate the performance of our method in detecting and estimating 
the location and velocity of multiple moving targets. 

We considered a scene of size 1100 x 1100 m 2 with flat topography centered at 
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Figure 7: Reconstructed reflectivity images (a) with the estimated velocity, v x = 
[—0.4, 6]m/s; and (b) with the correct target velocity, i.e., v x = [0, 6.2]m/s. The black 
circle indicates the true target position at time t = 0. 



Table 3: Analysis of the reflectivity images reconstructed using erroneous velocities 



Velocity error 
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tion shift (m) 
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No (Correct velocity) 








(65,65) th 


1 


|A< = 0.5m/s 


22 (|Az ± |) 


21.07 (|Az ± |) 


(76, 59) th 


0.6145 


|Av^| = 0.5m/s 


52 (|Az r |) 


42.14 (\Az r \) 


(62,39) th 


0.7237 



[ll,ll,0]km. The scene was discretized into 128 x 128 pixels, where [0,0, 0]m and 
[1100, 1100, 0]m correspond to the pixels (1,1) and (128,128), respectively. Fig. 10 
shows the scene with a static extended target and multiple moving targets along with 
their corresponding velocities. 

We assumed that the transmitter and receiver were traversing a circular trajectory 
given by 7c(s) = (11 + llcos(s),ll + 11 sin(s), 6.5) km. Let Jt(s) and 7#(s) 
denote the trajectories of the transmitter and receiver. We set 7t(s) = 7c(s) and 




Figure 8: (a) Reconstructed reflectivity image when Vh = [0, 6.7]m/s with a radial 
velocity error, Av£ = 0.5m/s. The black circle indicates the true target position at 
time t — 0. (b) X profiles, and (c) Y profiles. Solid lines show the X and Y profiles of 
the image reconstructed using the correct velocity shown in Fig. 7(b). 



1r(s) = 7c(s — §)• Note that the variable s in 7c is equal to j^t where V is the speed 
of the receiver or the transmitter, and R is the radius of the circular trajectory. We 
set the speed of the transmitter and receiver to 261 m/s. Fig. 11 shows the 3D view 
of the transmitter and receiver trajectories and the scene. 

The length of the signal was set to = 0.1707s. The circular trajectory was 
uniformly sampled into 2048 points, corresponding to f s = 7. 7339Hz. 

We assumed that the velocity of the targets is in the range of [—20, 20] x 
[— 20, 20]m/s and discretized the target velocity space into a 41 x 41 grid with the 
discretization step equal to lm/s. Thus, the velocity estimation precision is lm/s in 
our simulations. We reconstructed p Vh {x) images via the FBP method as described 
in Section 4.1. 




Figure 9: (a) Reconstructed reflectivity 
tangential velocity error, Av, 1 = 0.5m/s. 
position at time t = 0. (b) X profiles and 
profiles of Fig. 7(b). 



image when Vh = [0.5,6.2]m/s with a 
The black circle indicates the true target 
(c) Y profiles with a comparison with the 



We form the contrast image as described in Section 4.4. The result is shown 
in Fig. 12. We see from Fig. 12 that there are four dominant peaks marked 
with red circles. This indicates that there are four different velocities associated 
with the moving target scene. The velocities where the peaks are located are 
[—10, 15, 0], [0, 10, 0], [15, —5, 0]m/s and [0, 0, 0]m/s. The estimated velocities are equal 
to the true target velocities used in the simulations. 

Fig. 13 presents the reconstructed reflectivity images corresponding to the 
estimated velocities, i.e., [— 10, 15, 0], [0, 10, 0], [15, — 5, 0]m/s and [0,0, 0]m/s. We see 
that the targets are well-focused in the images formed using the correct velocity 
associated with each target. Note that Fig. 13(a) is the image reconstructed with 
v/j = [0,0, 0]m/s. In this case, the moving target imaging method described here is 
equivalent to the static target imaging method that we introduced in [27]. As expected 
only the static target is reconstructed in Fig. 13(a). 
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Figure 10: The moving target scene considered in the numerical simulations. 




Figure 11: A 3-D illustration of the simulation setup. The dark region denotes 
the scene (moving targets are not displayed) considered in the simulations. The 
transmitter/receiver antennas traverse a continuum of positions along the circular 
trajectory as shown by the red line. At a certain time instant, the transmitter and 
receiver are located at the positions indicated by the solid dots. 

8. Conclusion 

We have introduced a novel method for the synthetic aperture imaging of moving 
targets using ultra-narrowband transmitted waveforms. Starting from the first 
principle, we developed a novel forward model by correlating the received signal 
with a scaled version of the transmitted signal over a finite time interval. Unlike 
the conventional wideband SAR forward model, which is based on the start-stop 
approximation and high resolution delay measurements, this model does not use start- 
stop approximation and is based on the temporal Doppler induced by the movement of 
antennas and moving targets. The analysis of the forward model shows that the data 



Synthetic Aperture Imaging of Moving Targets using Ultra- Narrowband CW 



30 




-20 -15 -10 -5 5 10 15 20 



Figure 12: The contrast-image obtained from the reflectivity images reconstructed 
using a range of hypothesized velocities. 



used for reconstruction is the projections of the phase-space reflectivity onto the four- 
dimensional bistatic iso-Dopplcr manifolds. We next developed a FBP-type image 
reconstruction method to reconstruct the reflectivity of the scene and used a contrast 
optimization method to estimate the velocity of moving targets. The reflectivity 
reconstruction involves backprojecting the correlated signal onto the two-dimensional 
cross-sections of the four-dimensional iso-Doppler manifolds which we referred to as 
the position-space iso-Doppler contours for a range of hypothesized velocities. We 
showed that when the hypothesized velocity is equal to the true velocity of a scatterer, 
the singularity is reconstructed at the correct position and orientation. The PSF 
analysis shows that the visible singularities reconstructed are those that are at the 
intersection of position-space bistatic iso-Doppler curves and bistatic iso-Doppler-rate 
curves corresponding to the correct target velocity. We designed the filter so that 
the strength of the singularities are preserved at the correct velocity. The resulting 
filter depends not only on the antenna beam patterns, but also on the hypothesized 
velocity of targets. Using the image contrast optimization, we estimated the velocity 
of moving targets from a stack of reflectivity images. 

We have analyzed the resolution of the reconstructed reflectivity images and the 
velocity resolution available in the data by analyzing the PSF of the imaging operator 
and the temporal Doppler bandwidth of the correlated data. Our analysis shows that 
both the reflectivity and velocity resolutions are determined by the temporal duration 
and the carrier frequency of the transmitted waveforms. These findings are consistent 
with the Doppler ambiguity theory of CW waveforms. Additionally, our analysis has 
identified various other factors, such as the relative velocity between the antennas 
and the moving targets, the range of the antennas to the moving targets, the bistatic 
angle and the variation in the ground topography, etc. that affect the reflectivity and 
velocity resolution. 

We have analyzed the error in reconstructed reflectivity images due to error in 
target velocity. Our analysis leads to several important results in moving target 
imaging. First, it shows that the position error primarily depends on the component of 
the velocity error in the antenna look direction and the projection of the velocity error 
onto the planes perpendicular to the look direction and the trajectories of the antennas. 
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Figure 13: Reconstructed images corresponding to the velocities: (a) Vh = [0, 0, 0]m/s; 
(b) v h = [-10,15,0]m/s; (c) v h = [0, 10, 0]m/s; (d) v h = [15, -5, 0]m/s. 



Secondly, our analysis explains the artifacts expected due to moving targets when the 
image is reconstructed under a stationary scene assumption. Finally, it shows that 
the position error in the backprojected data is small when the error in the estimated 
velocity is small. Additionally, our error analysis method can be easily applied to 
understand and analyze the positioning errors due to errors in antenna positions. 

We presented extensive numerical simulations to verify our theoretical analysis 
and to illustrate the performance of our imaging method. 

We considered the bistatic scenario where the transmitting and receiving antennas 
are sufficiently far apart. The results for the monostatic case can be deduced by simply 
setting the two antenna trajectories to be equal. 
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Our moving target imaging method can be easily extended to incorporate imaging 
of airborne-targets and complex target motion models. While in the current paper, we 
assumed that the velocity of each target remains constant throughout the synthetic 
aperture, the forward model and the image formation method can be extended to 
include higher order kinetic parameters. We leave the investigation of this topic as a 
future research. 

Our imaging method can be implemented efficiently by using fast backprojection 
algorithms [37,38] or fast Fourier integral operator computation methods [29,31], and 
by utilizing parallel processing on graphics processing units [39] . 

Although our imaging scheme was developed in a deterministic setting, it is also 
applicable when the measurements are corrupted by additive white Gaussian noise [40] . 
When a priori information for the scene to be reconstructed is available and additive 
noise is colored, FBP-type inversion method presented in this paper can be extended 
as described in [41]. 

Finally, while our primarily interest is in radar imaging, our method is also 
applicable to other similar imaging problems such as those that may arise in acoustics. 
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Appendix A. 



Let 7T(R) = (7T(fl).7T(fl).7T(fl)) 7 7t(*) = (iT(R)^T(R)'iT(R)) T and 

{vi,V2,V3) T . We write 

7t(s) - (x + vs) • (7t(s) - v) 



\jt(s) - (x + vs 
+ (7r 0) - (X2 + «2s))(7t(s) - v 2 ) 
+ (7r0) - (a* + «3»))(7r(*) " v s)} ■ 
Note that x = (x, t/j(x)),x = {x\, IE2). 



[(7r( s ) - ( x i + «i s ))(7t( s ) ~ vi) 

(A.l) 



Thus, 



d~f T {s) - (x + vs) • (7t(s) - v) 



ds 
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(7r(s) - v) 



7t(s) - (x + vs)| 2 \jt(s) - (x + vs) 
x [(1t(s) - (x + vs)) • (7 T (s) - v)] 

|7t(s) - (x + vs)| 
+ (7t(s) - v 2 f + ( 7 | (s) - (x a + w 2 s)) 7 t(s) 
+ (7t(s) - «s) 2 + (7t0) - (^3 + «3«))7f(«)] 



-1 



|7t(s) 



(X + vs j 

1 



|7r(s) - (x + vs)| 



[(7 T (s)-(x + vs))-(7 T (s)-v)] 2 

[|7t(s) - v| 2 + (tt(s) - (x + vs)) • 7 T (s)] 
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|7 T (s) - (x + vs) 

+ (7t(s) - (x + vs)) -7t(s) 
1 



|7t(s) - (x + vs) 



|( 7t (s) - v) ± | 2 + ( 7T (s) - (x + vs)) • 7T (s) 



(A.2) 



Or,r (s,x) 

where 

(7r(s)-v) ± = (7 T (s)-v)-( 7t (s)^(x + vs))((7 T (s)^x + vs)) • ( 7t (s) - v))(A.3) 
denotes the projection of the relative velocity between the transmitter and the moving 
target 7r(s) — v onto the plane whose normal vector is along 7t(s) — (x + vs), 
aT,r(s,x) denotes the projection of the transmitter acceleration 7r(s) along 

7t(s) — (x + vs). We see that the summation of the two terms of (A.2) is the total 
radial acceleration of the transmitter evaluated at s with respect to the moving target 
located at x + vs on the ground at s. 
Using (A.2) and (34), we obtain 

|( 7T (s) - v)J 2 + ( 7T (s) - (x + vs)) • 7T (s) 



i ( \ fo 

/ d (s,x,v) = — 

CO 



|7r(s) - (x + vs) 
1 



|(7fl(s) - v) J 2 + (jr(s) - (x + vs)) • 7fi (s) 



(A.4) 



|7r(s) - (x + vs) | 
where similar to (A. 3), 

(7k(«) - v)_l = (7a(*)- v) - (7r(s)^-"(x~+ vs))(( 7 r(s)^-"(x~+ vs)) • ( 7r (s) - v))(A.5) 
denotes the projection of the relative velocity between the receiver and the moving 
target 7r(s) — v onto the plane whose normal vector is along 7j?(s) — (x + vs). 



Appendix B. 

Using (64), we have 

H„ (s, x) = 27rV x / d (s, x, v ) 

The first-order partial differential of (A.l) with respect to X\ is given by 
9(7t(s) - (x + vs)) • (7t(s) - v) 



dfd/dxi 
dfd/dx 2 



(B.l) 



dxi 



-I 



|7r(s)- (x + vs)| 2 

( 7 ^(s) - Xl ) - (ffi + j^i* + ^L.^ 3 )( 7 3 ( a ) - (^xg) + »**)) 

|7t(s) - (x + vs) 

x ((7t(s) - (x + vs)) • (tt(s) - v))] 
d 2 i> 



1 



|7r(s) - (x + vs) 



\1tK s ) - Vl) - (,"5 h — — w l s 



d 2 x\ 



"(7t( s ) - (V'^l)^) + V3S ^(~Q2^ Vl 4 



8x28x1 

a 2 v 



u 2 s)(7 T (s) - u 3 ) 



dx2dx 



■v 2 ) 



(B.2) 
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Note that v 3 = V x ^(x) ■ [v lt V2] = §^v 1 + §^v 2 . 

Similarly, the first-order partial differential of (A.l) with respect to z 2 can be 
expressed as follows: 

<9(7t(s) - (x + vs)) • (7t(s) - v) 



dx 2 



-I 



\1t(s) - (x + vs)| 2 

-( 7 2( S )_ X2 )_(^ 



dx\dx2 



vis + Q2^v 2 s)(^(s) - (ip(x 1 ,x 2 ) +v 3 s)) 



|7t(s) - (x + vs) 

x ((7t(s) - (x + vs)) • (7r(s) - v))] 



|7r(s) - (x + vs)| 



dx\dx 



-V\S ■ 



d 2 ip 
d 2 x 2 



f2s)(7r(s) - v 3 ) 



Wc define 



and 



Hence 



D = 



D 



"(7t( s ) - {tp{xi,x 2 ) + v 3 s))( 



I dip(x)/dxi 
f dip(x)/dx 2 



u u d' i x l u i + 95^7^2 



dx\dx 2 d 2 x 2 



v 2 ) 



(B.3) 
(B.4) 

(B.5) 



9(7r(g)-(x+v3))-(7T(^)-v) 

9(7r (s) - (x+vs) ) • (7T (s) - v) 

da; 2 



-[L> + D 2 s] 



(7t(s) - v)_ 



D 2 • ( 7T (s) - (x + vs)) 



(B.6) 



|7t(s) - (x + vs) 
where 

(7t(s)-v)± = (7t(s)-v)-(tt(s) - (x + vs))(7t(s) - (x + vs)) • (7t(s) - v) .(B.7) 
Thus, using (34), applying the derivation in (B.2), (B.3) and (B.6) to each 
component of d(jn(s) — (x + vs)) ■ (tr(s) — v)/dx, we obtain 



27t/q 
co 



(7r(s) - v )_ 



(7fl(s) - v oU 



[D + £> 2 s] 

|_|7r(sJ - (x + v s) |7 fl (s) - (x + v s) 
+D 2 ■ [( 7T (s) v s)) + ( lR (s) v s))]} . (B.8) 

Note that D 2 in (B.8) is given by (B.5) with vi,v 2 replaced with no. 1,1*0.2 where 

Vo = [V0,1,V0,2, Va-V^) ' KbfO^]]. 
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Appendix C. 
Using (34), we have 

fd(s,z,v z +eAv z ) = — (7t(s) - (z + v z .s + eAv z s)) • ( 7 t(s) v z eAv z ) 

+(7fl(«) - (z + vTs + eAv z s)) • (<y fl («) v z eAv z ) . (C.l) 

Now we calculate the first derivative of (C.l) with respect to e. Let us consider 
the derivative of the first item in the square bracket in (C.l). 

d 



de 



(7t(s) - (z + v z s + eAv z s)) ■ (tt(s) 



eAv z ) 



= (7t(s) - (z + v z s + eAv z s))'| e=0 • (7t(s) v z eAv z ) 
+ (7t(s) - (z+^s + eAv z s)) ■ (-Av z )| e=0 
where []' denotes the first derivative with respect to e, 

(7t(s) - (z + v z s + eAv z s))'| £=0 
_ -[Avs - [( 7 t(s) ;r Tz + v.a)) • Av z s]( 7T (s) ^"(7+ v,*))] 



(C.2) 



|7t(s) - (z + v z s) 



-sAvd 



| 7 t(s) - (z + v z s)| 

Note that Av^ T is the projection of Av z on the plane whose normal direction is 
along(7 T (s) - (z + v z s)). 

Thus, substituting (C.3) into (C.2), and the result back into (C.2), we have 

d e fd(s,z,v z +eAt> z )| e=0 
fo 



S AvJ- r -(7r(s)-v z ) 

\1t(s) - (z + v z s)| 
sA^-R • ( 7fl (s) - v.) 



- Av z • ( 7 t(s) - (z + v z s)) 

- Av z • ( 7fi (s) - (z + v z s)) 



|7ii(s) - (z + v z s) 
We assume flat topography. From (69) of the manuscript, we have 



(C.4) 



V z f d (s,z,v z ) 



CQ 



D 



(7t(s)-v z )- 



(7fl(«) - v.)- 



|7t(s) - (z + v z s) | 7 _r(s) - (z + v z s) 



fo f 27t/ „ , \ 
-a v ^(s,z) 



C'O 



(C.5) 



Plugging (C.4) and (C.5) into (91), we obtain 

■Av^-OM-O-v.) , Av^.fe( S )-v z ) 



f.S 



|7t(s) - (z + v z s)| |7fl(s) - (z + v z s)| 
eAv z • [7t(s) - (z + v z s) + 7r(s) - (z + v z s)] 
Az ■ S Vz (s,z) . 



(C6) 
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Appendix D. 
Using (56), we have 



fd(s,z,v z + eAv z ) = — 



Let 



|(7t(s)-(v z + 6Av z )U| 2 
|7t(s)- (z + (v z + eAv z )s)| 

+(7t(s) - (z + K + eAv z )s)) • j T (s) 

| |( 7fi OO-(v z + e Av z )K| 2 

\-y R (s) - (z + (v z + eAv z )s)| 
+(7fl(«) - + K + eAv z )s)) ■ 7R (s) 

-|(7T( S )-(v + 6At;)K| 2 . 



|7t(s) - (z + (v + eAv)s)| 
Calculating the first derivative of (D.2) with respect to e, we obtain 



de 



1 



|7t(s) - (z + (v z + eAv z )s 
1 



|7t(s) - (z + (v z + eAv z )s 
where []' denotes the first derivative with respect to e 



|7 T (s)-(z + (v z +eAv z )s) 



_ 7tQ) ~ (z + v z g + eAv z s) ■ Av z s 
\jt{s) - (z + v z s + eAv z s)| 2 



and 



(|(7t(s) - (v z + £ Av z )) ± | 2 )' = 2|( 7T ( S ) - (v z + eAv.))±| 



x l(7T(s)-(v z +eAv z )) ± |'. 



In (D.5), 

|(7T(.s)-(v z +eAv z )) ± |' = 
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(D.l) 
(D.2) 



|(7T(s)-(v z + eAv z )) x p 



(|(7t( S )-(z + (v z + £ Av z )K| 2 )', (D.3) 



(D.4) 



(D.5) 



|(7r(s)-(v z + eAv z )) x | 
x (7t(s) - (v z + eAv z ))i • {j T (s) - (v z + eAv z )) ± 



(D.6) 



where 



(7t(s) - (v z + eAv z ))'_ L | e=0 



= -Av z - 



(7r(s) - (z + v z s)) • (tt(s) - v z ) 



-Av; 



_,T 



|7t(s) - (z + v z s)| 

7t (s) - (z + v z s) — -5- -• 

L|7t(s) - (z + v z s) 

Av ■ 7t(s) — (z + v,s) | 

s 7t(s) - (z + v z s)) • (tt(s) - v z ) 



(7t(s) - v z ) 



7t(s) - (z + v z s) 



|7t(s) - (z + v z s)| 
|7t(s) - (z + v z s)| 



(D.7) 
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where Av^ ,T,fl is the projection of Av z on the plane whose normal direction is along 

(1tm{s) - (z + v„s)). 

Substituting (D.7) into (D.6) and then the result back into (D.5), we have 

(|(7T(s)-(v + eAu))J 2 )'| e=0 

S7t(s) - (z + v z s)) • (7t(s) - v z ) 



= -2Ai 



(7r(s) - v z ) ± 



1 - 



|7t(s) - (z + v z s) 



(D.8) 



Using (D.8), (D.4) and (D.3), we obtain 



— e=0 = Av z -7t s - (z + v z s — — — -2 

oe 7t(s) - (z + v z s) 2 



Av^ • (j T (s) - v.)x 



S7t(s) - (z + v z s)) • (tt(s) - v z ) 



|7t(«) - (z + v - zs) 
Thus, using (D.9) and (C.3), we have 
dj d (s, z,v + eAv)\ e=0 



fo 

CO 



s7t(s) 



|7r(s) - (z + v z s) 



2(7t(s) - v z )_l 



(D.9) 



|7t(s) - (z + v z .s) |7 T (.s) - (z + v z s; 
s|(7r(s) - v z ) ± | 2 



-C T (z, v z , s) 



-Av z • 7t(s) - (z + v z s) 



-A 



S7fl(s] 



z + v z sjp 
2(7«(a) - v.) 



|7r(s 

, 2( 7fl ( s ; 

|7ft(«) - 
s l(7/?( s ) - v z)±| 2 



l7fl(«) - (z + v z .s) hB.(s) - (z + v z s 
Av z -7fl(s) - (z + v z s) 



^C fl (z,v z ,s) 



where 



\jr(s) - (z + v z s)| 2 _ 
s Jt,r( s ) - (z + v z s) • (7r,fl(s) - v z ) 



C T ,ii(z,v z ,s) = 1 - 

|7t,h(s) - (z + v z s)[ 

Now let us consider V z fd(s, z,v z ). We assume flat topography. We have 
Vz ( \(7t(s)-v z ) ± \ 2 
Z \ \Jt{s) - (z + v z s)| 
I(7t( S )-v z ) ± | 2 



(D.10) 



(D.ll) 



V 2 



|7t(s) - (z + v z ,s 
1 



|7r(s) - (z + v z s 
1 



+ 7t(s) - (z + v z s) • 7r(s) 

+ V z (7t(s) - (z + v z s) • 7t(s)) 
|(7t(s) - v z ) ± | 2 



V z |(7t(s) - v z )_l| 



|Tt(s) - (z + v z s 
+(V z 7t(s) - (z + v z s)) • 7t(s) 



In (D.12), 
V 



1 _ D • 1t(s) - (z + v z s) 

|7t(«) ~ (z + v z s) |7t(«) - (z + v z s)| 2 



(D.12) 
(D.13) 
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7t(s) - (z + v z s) 



V z 7t(s) - (z + vs) 



|7t(s) - (z + v z s)| 
1 



|7t(s) - (z + vs) 



[D • 7 t(s) 



where for flat topography, 

1 



L> = 



1 



and 



where 



V*|(7t(s)-v z ) x 

V z (7t(s) - v z )_l 

7t(s) - (z + v z s) 
|7t(s) - (z + v z s) 
(7t(s) - v z ) • 7 T (s) - (z + vs) 



2V z (7t(s) - v z )_l • (7t(s) - v z )_l 



(D.14) 

(D.15) 
(D.16) 



|7t(s) - (z + v z s)| 
7t(s) - (z + v z ,s 



[D-jt(s) - (z + v z s)](tt(s) - v z ) -7t(s) - (z + v z s) 



|7t(s) - (z + v z s) 
7t(s) - (z + v z s) 



|7r(s) - (z + v z s 
Substituting (D.17) into (D.16), we obtain 



[D ■ 7 T (s) - (z + v z s)](7 T (s) - v z ) • 7 T (s) - (z + v z s) 
( j D-7t(s)-v z ). (D.17) 



V z |( 7t (s)-v)J 2 = 2 



n (tt (s) - v z ) -7t(s) - (z + vs) 
|7t(s) - (z + vs) 



D-(7t(s)-vK (D.18) 



Using (D.18), (D.13) and (D.14), (D.12) becomes 
^ ( |(7t( S )-v z ) ± | 2 



V|7r(s) - (z + v z s 

2 (7t(s) - v z ) • 7 T (s) - (z + v z s) 



|7x(s) - (z + v z s)| 2 
I(7t( S )-v z ) ± | 2 



+ 7t(s) - (z + v z s) • 7t(s) 

(7t(s) - v z )_l 



| 7t (s) - (z + v z s)| 2 

j£(g) 

|7t(s) - (z + v z s)| 



7t(s) - (z + v z s) 



(D.19) 



Using (D.19), (D.10), considering the fact that V z (d s f d ) = d s (\7 z f d ) = ^-d s S Vx , 



after rearrangement, we obtain 



S7t(s) 



2(7t(s) - v z 



|7t(s) - (z + v z s)| |7 T (s) - (z + v z s; 

|(7t(s)-v z )J 2 



+esAv z • 7t(s) - (z + v z s 



-eAv 



|7t(s) - (z + v z s)| 2 
sJr(s) 2(7r(s)-v z )x 



|7r(«) - (z + v z s)| |7i?( s ) ~ ( z + v zs) 



f C T (z,v z ,s) 



Cr(z,V z ,s) 
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^Av z . 7 , (S )^r + v zS) i(-fr«-v.)xi a 



|7r(s) - (z + v z s)| 2 



-Az-S^s,*)-^ (D.20) 



2tt/ ( 



o 



where 



2^/o U7t(s)-(z + v z s)| 2 

|(7k(*)-v.)±| 2 




|7b(s) - (z + v z .s)| 2 

2 (7r(s) - v z ) • 7r(s) - (z + v z s) 
|7t(s) - (z + v z s)| 2 

(7g(g) - v z ) • 7fl(£) - ( z + v zs) , . , s v 

+ 2 i ri — ? — I V2 (7Ms) - v z )_ 

|7h(s) - (z + v z s)| 2 

7t( s ) 7a(s) 



|7t(s) - (z + v z s)| |7h(s) - (z + v z s)| 



(D.21) 
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